Large fluctuations in stochastic population dynamics: momentum space calculations 
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Abstract. Momentum-space representation renders an interesting perspective to the- 
ory of large fluctuations in populations undergoing Markovian stochastic gain-loss processes. 
This representation is obtained when the master equation for the probability distribution 
of the population size is transformed into an evolution equation for the probability gen- 
erating function. Spectral decomposition then brings about an eigenvalue problem for a 
non-Hermitian linear differential operator. The ground-state eigenmode encodes the sta- 
tionary distribution of the population size. For long-lived metastable populations which 
exhibit extinction or escape to another metastable state, the quasi-stationary distribution 
and the mean time to extinction or escape are encoded by the eigenmode and eigenvalue of 
the lowest excited state. If the average population size in the stationary or quasi-stationary 
state is large, the corresponding eigenvalue problem can be solved via WKB approximation 
amended by other asymptotic methods. We illustrate these ideas in several model examples. 
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I. INTRODUCTION 

This work deals with dynamics of populations experiencing intrinsic noise caused by the discreteness of individuals 
and stochastic character of their interactions. When the average size N of such a (stationary or quasi-stationary) 
population is large, the noise-induced fluctuations in the observed number of individuals are typically small, and only 
rarely large. In many applications, however, the rare large fluctuations can be very important. This is certainly true 
when their consequences are catastrophic, such as in the case of extinction of an isolated self-regulating population 
after having maintained a long-lived metastable state, with applications ranging from population biology [1, 2] and 
epidemiology [1, 3] to genetic regulatory networks in living cells [4]. Another example of a catastrophic transition driven 
by a rare large intrinsic fluctuation is population explosion [5]. Rare large fluctuations may also induce stochastic 
switches between different metastable states [6]; these appear in genetic regulatory networks [7] and in other contexts. 
Less dramatic but still important examples involve large fluctuations in the production rates of molecules on the 
surfaces of micron-sized dust grains in the interstellar medium, where the number of atoms, participating in the 
chemical reactions, can be relatively small [8, 9]. As stochastic population dynamics is far from equilibrium and 
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therefore defies standard methods of equilibrium statistical mechanics, large fluctuations of stochastic populations are 
of much interest to physics [10, 11]. 

In this paper we consider single-species populations which are well mixed, so that spatial degrees of freedom are 
irrelevant. To account for stochasticity of gain- loss processes (in the following - reactions) and discreteness of the 
individuals (in the following - particles) , we assume a Markov process and employ a master equation which describes 
the time evolution of the probability Pn{t) to have a population size n at (continuous) time t. If the population 
exhibits neither extinction, nor a switch to another metastable state or to an infinite population size, a natural goal 
is to determine the stationary probability distribution of the population size [12]. For metastable populations - 
populations experiencing either extinction, or switches between different metastable states - one is usually interested 
in the mean time to extinction or escape (MTE), and in the long-lived gwasi-stationary distribution (QSD) of the 
population size. For single-step processes these quantities can be calculated by standard methods. The MTE can be 
calculated exactly by employing the backward master equation [10, 11]. This procedure yields an exact but unwieldy 
analytic expression for the MTE which, for a large population size in the metastable state, can be simplified via a 
saddle-point approximation [13]. In its turn, the QSD of a single-step process can be found, in many cases, through 
a recursion. 

For multi-step processes neither the MTE, nor QSD can be calculated exactly. Many practitioners have used, in 
different contexts of physics, chemistry, population biology, epidemiology, cell biology, etc, what is often eallcid "the 
diffusion approximation" : an approximation of the master equation by a Fokker-Planck equation. The latter can be 
obtained via the van Kampen system size expansion, or other related prescriptions. With a Fokker-Planck equation 
at hand, the MTE and QSD can again be evaluated by standard methods [10, 11]. Unfortunately, this approximation 
is in general uncontrolled, and fails in its description of the tails of the QSD. As a result, it gives exponentially large 
errors in the MTE [13-16]. 

Until recently, the MTE and QSD had been calculated accmately only for a few model problems involving multi- 
step processes. Recently, Escudero and Kamenev [17] and Assaf and Meerson [18] addressed quite a general set of 
reactions and developed controlled WKB approximations for the MTE and QSD for population switches [17] and 
population extinction [18]. When necessary, the WKB approximation must be supplemented by a recursive solution 
of the master equation at small population sizes [5, 18], and by the van Kampen system size expansion in narrow 
regions where the "fast" and "slow" WKB modes are coupled and the WKB approximation fails [5, 17, 18]. 

The techniques developed in Refs. [17] and [18] (see also Refs. [5, 6, 15]) were formulated directly in the space of 
population size n. An alternative approach invokes a complementary space which can be interpreted as a momentum 
space. The momentum-space representation is obtained when the master equation - an infinite set of linear ordinary 
differential equations - is transformed into a single evolution equation - a linear partial differential equation - for 
the probability generating function G{p,t). Here p, a complementary variable, is conjugate to the population size n 
and plays the role of "momentum" in an effective Hamiltonian system which encodes, in the leading order of 1/7V- 
expansion, the stochastic population dynamics. One can then perform spectral decomposition of this linear partial 
differential equation for G{p,t). In order to describe the stationary or metastable states, it suffices to consider the 
ground state and the lowest excited state of this spectral decomposition, whereas higher modes only contribute to 
short-time transients [16, 19]. 

The ordinary differential equations for the ground state and the lowest excited state are determined by the specific 
set of reactions the population undergoes. The order of these equations is equal to the highest order of inter-particle 
reactions. For example, for two- (three-) body reactions the equations are of the second (third) order, etc. In general, 
these ordinary differential equations cannot be solved exactly, and some perturbation techniques, employing the small 
parameter 1/iV <C 1, need to be used. 

The momentum-space spectral theory was developed [16, 19, 20] for two-body reactions. Here we extend the theory 
to any many-body reactions. We also determine, in the general case, the previously unknown boundary conditions 
for the above-mentioned eigenvalue problems. If there is no absorbing state at infinity, the boundary conditions are 
"self-generated" by the demand that the probability generating function G{p, t) be, at any t, an entire function on 
the complex plane p [21]. We show that, for two-body reactions, the population extinction problem can always be 
solved by matching the exact solution of a quasi-stationary equation for the lowest excited state (see below) with 
a perturbative solution of a non-quasi-stationary equation for the same state. This procedure always works when 
N is sufficiently large. For three-, four-, . . . body reactions the spectral decomposition can be used in conjunction 
with a p-space WKB (Wentzel-Kramers-Brillouin) approximation which employs the same small parameter 1/N but 
does not rely on exact solution of the quasi-stationary equation. We find that there is a region of p where the WKB 
approximation breaks down, and a region of p where its accuracy is insufficient. In the former region a boundary-layer 
solution can be found and matched with the WKB solution. In the latter region a simple non- WKB perturbative 
solution can be obtained. The theory extensions presented here turn the momentum-space spectral theory of large 
fluctuations into a more general tool. 

As the evolution equation for G{p, t) is equivalent to the master equation, the p-space approach is clearly advanta- 
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gcous. compared to the n-space approach, when the problem in the p space admits an exact solution, see Refs. [8, 11]. 
Otherwise, the technical advantages of the p-space approach are not a priori obvious. In any case, it provides a viable 
alternative, and an interesting perspective, to theory of large fluctuations of stochastic populations. 

Here is the layout of the rest of the paper. Section II briefly introduces the momentum-space spectral formalism, 
whereas in sections III and IV we describe the methods of solution and illustrate them on several model examples. 
Sec. Ill deals with a well-studied prototypical chemical reaction scheme which describes a stationary production of 
hydrogen molecules on interstellar dust grains. Here we show that the WKB approximation not only gives accurate 
results for the production rate (including its fluctuations) of hydrogen molecules, but also yields a complete stationary 
probability distribution function of the number of hydrogen atoms, including its non-Gaussian tails. In Sec. IV we 
deal with isolated populations undergoing intrinsic-noise-driven extinction after maintaining a long-lived metastable 
state. Here, after some general arguments, we consider two diff'erent examples - one studied previously and one new - 
and determine the MTE and QSD. Throughout the paper we compare our analytical results with numerical solutions 
of the pertinent master equation and, when possible, with previous analytical results. Section V summarizes our 
findings and discusses the advantages and disadvantages of the p-space method compared with the "real" space WKB 
method [5, 6, 15, 17, 18]. 



II. MASTER EQUATION, PROBABILITY GENERATING FUNCTION AND SPECTRAL 

FORMULATION 



Populations consist of discrete "particles" undergoing stochastic gain and loss reactions. To account for both 
discreteness and stochasticity, we assume the Markov property, see e.g. Refs. [10, 11], and employ the master 
equation 

Pn{t) = Wn'nPw ' Wnn' Pn (1) 

which describes the time evolution of the probability distribution function P„(t) to have n particles at time t. Here 
Wnn' is the transition rate matrix; it is assumed that Pn<o = 0. 
The probability generating function, see e.g. Refs. [10, 11], is deflned as 

oo 

G{p,t) = J2p"Pn{t). (2) 
n=0 

Here p is an auxiliary variable which is conjugate to the number of particles n. Once G{p, t) is known, the probability 
distribution function Pn{t) is given by the Taylor coefficients 



put) 



nl dp" 

or, alternatively, by employing the Cauchy theorem 



(3) 

p=0 



where the integration has to be performed over a closed contour in the complex p-plane around the singular point 
p = 0. For stochastic populations which do not exhibit population explosion [5], the probability P„(t) decays faster 
than exponentially at large n. Therefore, G{p,t) is an entire function of p on the complex p-plane [21]. 

If the reaction rates are polynomial in n, one can transform the master equation (1) into a single linear partial 
difli'erential equation for the probability generating function, 

§=CO. (5) 

where £ is a linear differential operator which includes powers of the partial differentiation operator d/dp. Equation (5) 
is exact and equivalent to the master equation (1). If only one-body reactions are present, C is of flrst order in d/dp, 
and Eq. (5) can be solved by characteristics [11]. For many-body reactions one can proceed by expanding G{p,t) in 
the yet unknown eigenmodes and eigenvalues of the problem [16, 19, 20]: 

oo 

G{p, t) = G^tip) + ^ afe</.fe(p)e-^'=* . (6) 
fe=i 
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As a result, partial differential equation (5) is transformed into an infinite set of ordinary differential equations: for the 
(stationary) ground-state mode Ggtip) and for the cigcnmodcs of excited states {4'kip)}'kLi- By virtue of Eq. (3) or 
(4), the ground state eigenmode determines the stationary probability distribution function of the system. If a long- 
lived population ultimately goes extinct, the stationary distribution is trivial: P„ = Sn,o, where dn,o is the Kronecker's 
delta. What is of interest in this case is the quasi-stationary distribution and its (exponentially long) decay time which 
yields an accurate approximation to the MTE. These quantities are determinc;d by the lowest excited eigenmode 4'{p) 
and the eigenvalue Ei, respectively [16, 19]. The higher modes only contribute to short-time transients. Therefore, 
in the following we will focus on determining Gst or solving th e eigenvalue problem for (p) and Ei . 



III. STATIONARY DISTRIBUTIONS: GROUND-STATE CALCULATIONS 

As a first example, we consider a simple model of production of H2 molecules on micron-sized dust grains in 
interstellar medium. This model was investigated by Green et. al. [8], who computed the stationary probability 
distribution function of the number of hydrogen atoms via finding an exact solution to the ordinary differential 
equation for Gst- The same results were obtained, by a different method, by Biham and Lipshtat [9]. We will use 
this problem as a benchmark of the ground-state calculations using the momentum-space WKB approach. As we will 
see, this approach gives, for A'' 1, an accurate approximate solution for Gst{p), and so it can be employed for many 
other models where no exact solutions are available. 

Consider the following set of reactions: absorption of il-atoms by the grain surface $ H, desorption of H- 

atoms, 7? — )• 0, and formation of il2-niolecules from pairs of H-atovas which can be formally described as annihilation 
4 0. 

To calculate the production rate of il2-molecules, one needs to determine the stationary probability distribution 

fimction of the iJ-atoms, P„(i — > 00). For convenience, we rescale time and reaction rates by the desorption rate /3 and 
denote N = 2/3/7 and R = a'y/{2(3'^). Ignoring fluctuations, one can write down the following (rescaled) deterministic 
rate equation: 

2 

n = NR-n- —fi^ , (7) 

where n{t) ':$> 1 is the average population size. The only positive fixed point of this equation, 

N 



n = —{VT+8R-1), (8) 

is attracting, and the stationary probability distribution function P„ is expected to be peaked around it. The master 
equation describing the stochastic dynamics of this system in rescaled time is 

j^Pn{t) = ^[{n + 2)(n + l)Pn+2{t) - n{n - l)P„(t)] + [{n + l)P„+i(t) - nP„(t)] + AfP(P„_i - P„) . (9) 
This yields the following partial differential equation for G{p, t) [8] : 

The steady-state solution G^t obeys the ordinary differential equation 

^ (1 + p)G':, + - NRGst = , (11) 

where primes denote the p-derivatives. The boundary conditions are "self-generated" . Indeed, equality G{p = l,t) = 1 
holds at all times. This reflects conservation of probability, see Eq. (2). Therefore, 

G,t(l) = l. (12) 

Furthermore, Eq. (11) has a singular point at p = —1. As Gst{p) must be analytic at p = —1, we demand 

G;,(-l)-AfPG,i(-l) = 0. (13) 

The boundary-value problem (11)-(13) is exactly solvable in special functions [8]. For a general set of reactions, 
however, one cannot expect an exact solution. Still, one can employ the small parameter 1/N to develop an accurate 
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analytical approximation. To illustrate this point we will proceed as we were unaware of the exact solution, and then 
compare the approximate solution with the exact one. As the small parameter 1/N appears in the coefficient of the 
highest derivative, it is natural to use (a dissipative variant of) the stationary WKB approximation in the f)-space 
[20]. The WKB ansatz is 

G.t(p) = a(p)e-^^W, (14) 

where the action S{p) and amplitude a{p) are non-negative functions of p. Using this ansatz in Eq. (10) with a zero 
left hand side, we obtain 



1 

TV 



(1 - p^) [o" - 2NS'a' - NS"a + N'^{S'fa] + (1 - p){a' - NS'a -NRa) = 0. (15) 



In the leading order 0{N) we obtain a stationary Hamilton- Jacobi equation H\p, —S'{p)] = with zero energy, c/. 
Ref. [22]. The effective Hamiltonian is 

H{p,q) = {l-p)[{l+p)q^ + q-R], (16) 

where we have introduced q{p) = —S'{p): the reaction coordinate conjugate to the momentum p. The trivial zero- 
energy phase orbit p = 1 is an invariant line of the Hamiltonian; it corresponds to the deterministic dynamics [22]. 
Indeed, the Hamilton's equation for q, 



q = R-q-2q^, 

coincides, in view of the relation q = n/N , with the deterministic rate equation (7). Hamiltonian (16) also has two 
nontrivial invariant zero-energy lines which are composed of the two solutions, q-{p) and q+{p), of the quadratic 
equation (1 + p)q'^ + q — R = 0: 

I \ -l-^^(P) / N -1+U(P) 

^-(^) = ^(rTi) ' ^+(^^ = ^(rTi) • ^^'^ 

Here we have denoted 



v(3>) = ^/l + 4R{l+p). (18) 

The phase plane of this system is shown in Fig. 1. The phase orbits q = q-{p) must be discarded. This is because 

q-{p) diverges at p = —1, whereas Gst{p), and therefore S{p), must be analytic everywhere. 

The remaining nontrivial zero-energy phase orbit q+{p) = q{p) has a special role. It describes the most probable 
path along which the system evolves, (almost) with certainty, in the course of a fluctuation bringing the system from 
the fixed point {l,qi) in the phase space {p,q) to a given point, see Fig. 1. Here qi = (l/4)(\/l -|- 8-R — 1)] is the 
attracting point of the deterministic rate equation, see Eq. (8). 

Integrating the equation S'{p) = — g+(p), we obtain 

5(p) = -^(p) + z;(l)+ln^M±l, (19) 

where we have fixed the definitions of a{p) and S{p) by demanding S{p = 1) = 0. 

To calculate the amplitude a{p) we proceed to the subleading 0(1) order in Eq. (15): 

-2(1 + p)S'a' - (1 + p)S"a + a' = 0. (20) 

Using S{p) from Eq. (19), we arrive at a first-order ordinary differential equation for a{p), 

a'{p) _ 4i?2(i+p) 



a{p) t;(p)2[l -h u(p)]2 
Solving this equation, we obtain the WKB solution 



(21) 



WKB^^^ _ ^(1)^/^ [^ + V(P)] -NS(P) 



v{pY/^[l + v{l)f ' ^^^^ 

where the integration constant is chosen so as to obey boundary condition (12). As one can easily check, WKB 
solution (22) also obeys boundary condition (13). 
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FIG. 1: Molecular hydrogen production on a grain. Shown are zero-energy orbits of Hamiltonian (16) on the phase plane (p, q). 
The thick solid line corresponds to the instanton q — q+{p), see Eq. (17). The motion along the vertical line p = 1 is described 
by the deterministic rate equation (7). The dashed lines depict the branch q — q-{p). It is non-physical at g < and does not 
contribute to the WKB solution at g > 0. 



As expected, pre-exponent a{p) of the WKB solution (22) diverges at the turning point ptp = — 1 — l/(4i?) < —1 of 
the zero-energy phase orbit, see Fig. 1. As a result, the WKB solution breaks down in a close vicinity of this point. 
At p < ptp a WKB solution of a different nature appears: it exhibits decaying oscillations as a function of p. The 
oscillating WKB solution can be found by treating S{p) as a complex-valued, rather than real, function. We will not 
need the oscillating solution, because the non-oscillating one, Eq. (22), turns out to be sufficient for the purpose of 
calculating the probabilities Pm see below. 

Now we can compare WKB solution (22) with the exact solution of the problem (11)-(13), derived by Green et al. 
[8]: 



^' \l+p) /jv-i(27VV2i?) 



where Ik{w) is the modified Bessel function. To this end let us calculate the large- asymptote of /Ar-i[2A^-\/i?(l + p)] 
by using the integral definition of the modified Bessel function [23] 



where r(. . .) is the Euler Gamma function. As A^ 3> 1, we can evaluate the integral by the saddle point approximation 
[24]. Denoting /(<) = ln(l — i^) — 2-y/i?(l + p)t, we find the relevant saddle point 



l-^l + ARjl+p) _ _ I vjp) - I 
2/R(TT^ \v{p) + l' 

with v{p) from Eq. (18). Then, expanding f{t) ~ f{t^) + (1/2) f"{t^){t - t^f with /"(t,) = -v{p)[l + v{p)], and 
performing the Gaussian integration, we obtain the N ^ 1 asymptote 

In^i[2N^R{1+p)] ~ — ^^""^^^ , [Ar2i?(i +p)]^ e^{^(p)-i+'"2-'"[i+"(f)l> . (25) 

2V2r(Ar-l/2)v/A^ ^ V 

Note that the saddle point approximation is valid on the entire segment ~1 < P < 1- In particular, Eq. (25) with 
p — \ yields the A^ ^ 1 asymptote of the denominator of Eq. (23). Now one can see that the large- A^ asymptote 
of Eq. (23) exactly coincides with WKB solution (22). Actually, the WKB result is indistinguishable from the exact 
result already for A^ = 10, see Fig. 2. 
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FIG. 2: Molecular hydrogen production on a grain. Shown is a comparison of WKB result (22) for Gst(p) (dashed line) and 
exact result (23) (solid line) for A*' = 10 and R= 1. The agreement is excellent even for this moderate A''. 

Of a primary interest in the context of astrochemistry is the mean and variance of the steady-state production rate 
of H2 molecules. Going back to physical units, we can write the mean steady-state production rate as 

00 

n=0 



n{H2) 



W) + 1] 



1 - 



7Vt;2(l) 



(26) 



where v{jp) is given by Eq. (18). One can check that this expression coincides with that obtained from the exact 
result, see Eq. (22) in Ref. [8], in the leading- and subleading-order at ^ 1. The leading term in Eq. (26) is what 
the deterministic rate equation (7) predicts. 

Now consider the variance of the steady-state production rate of H2 molecules: 



Using identity 

n^{n- if 
we obtain the exact relation 



V{H2)^^[{n\n-lf)-{n{n-l))^] 



n{n - l){n - 2){n -3)+ An{n - l){n - 2) + 2n{n - 1) , 



V{H2) = I {G^^(l) + 4G"'(1) + 2G"(1) - [G"(l)]2} 

From WKB solution (22) we obtain in the leading order 

UjN^R^ [v{l) + 6R+l] 



V{H2) 



v{i) b(i)+ir 



(27) 



The relative fluctuations of the production rate, \/V /TZ, scale with N as N^^^^, as expected. 

Actually, the WKB approximation yields the whole stationary probability distribution function of the number of 
H atoms. Green et. at [8] obtained this distribution exactly from Eqs. (23) and (3): 



„ _ ^^{N^RT/^ lN+n-l{2NVR) 
— ^ ^ i ; 



{2NV2R) 



(28) 
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FIG. 3: Molecular hydrogen production on a grain. Shown is the natural logarithm of the stationary distribution P„ versus n 
for TV = 50 and R = 1. The solid line is WKB approximation (29), the dashed line is exact solution (28), and the dash-dotted 
line is the Gaussian approximation. The WKB approximation and the exact solution are indistinguishable for all n. The 
non-Gaussian tails of the distribution cannot be described correctly by the van Kampen system size expansion. The inset 
shows, by different symbols, the small-n asymptote of the distribution obtained analytically and numerically. 



The 3> 1, 71 S> 1 asymptote of (28) can be written as 



Pn 



V(T+9RT) 1 + uiq) 



.Af{ln[l+t,(l)]-i,(l)+9-H(l+<?)«(g)-ln[(l+g)(l+«(g))]-gln[<,(l+g)(l+«(g))/(2fl)]} 



v/27rgiVu(g) l + v{l) 
where q = n/N, v{p) is given by Eq. (18) and 



u{q) = v/1 + 4i?/(l + g)2 . 
Now we compare Eq. (29) with the WKB result, obtained from Eqs. (4) and (22): 



jWKB 



dp 



i;(l)i/2 [l + v{p)] 



2m J ^pw(p)i/2[i + i;(i)] 



'N S(p)-n\np 



(29) 



(30) 



(31) 



where S{p) is given by Eq. (19). As n 3> 1, we can evaluate the integral via the saddle point approximation. Let 
f{p) = —NS{p) — nlnp. The saddle point is at = g(l + q)[l + u{q)]/{2R), where u{q) is given by Eq. (30). As 
f"{p*) > 0, the integration contour in the vicinity of the saddle point must be chosen perpendicular to the real axis. 
This adds an additional phase of e*'^/^ to the solution [24], which cancels i in the denominator of Eq. (31). After the 
Gaussian integration and some algebra Eq. (31) coincides with Eq. (29). Finally, one can calculate P„ at iV ^ 1 but 
n = 0(1) by directly differentiating the WKB result (22) for Gst{p), see Eq. (3). The resulting probability distribution 
function is shown in Fig. 3. As one can see, the agreement between the WKB distribution and the exact distribution 
is excellent for all n. 



IV. METASTABILITY AND EXTINCTION: FIRST-EXCITED-STATE CALCULATIONS 

Now we switch to isolated stochastic populations, so that there is no influx of particles into the system. If there is 
no population explosion, isolated populations ultimately undergo extinction with probability one. The deterministic 
rate equation for such a population can be written as 

n = nil^{n), (32) 

where ipln) is a smooth function. In the following we assume ip{0) > 0, so that n = is a repelling fixed point of 
Eq. (32). The deterministically stable population size corresponds to an attracting fixed point n — ni > 0. According 
to the classification of Ref. [18], such populations exhibit scenario A of extinction. 
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Let Til = 0{N) ^ 1. After a short relaxation time tr, the population typically converges into a long-lived metastable 
state whose population size distribution is peaked around n = ni. This metastable probability distribution function 
is encoded in the lowest excited eigenmode (j){p) = 4>i{p) of the probability generating function G{p,t) (6). Indeed, at 
t » tr, the higher eigenmodes in the spectral expansion (6) have already decayed, and G{p,t) can be approximated 
as [16, 19] 

G(p,t)~l-0(p)e-^*, (33) 

where the lowest excited eigenfunction is normalized so that ^(0) = 1. The (exponentially small) lowest excited eigen- 
value E = El determines the MTE of the population, E ~ t~J' . The slowly time-dependent probability distribution 
function of the population size, at t ^ tr, is 

P„>o(t) ^ 7r„e-*/-- , Po(i) - 1 - e"*/"'^ • (34) 

That is, the metastable probability distribution function exponentially slowly decays in time, whereas the extinction 
probability Po{t) exponentially slowly grows and reaches 1 at t — > oo. The shape function 7r„ of the metastable 
distribution is called the quasi-stationary distribution (QSD). The QSD and MTE of a metastable population can be 
obtained by solving the eigenvalue problem for (l){p) and E, respectively. We now discuss some general properties of 
the solution to this eigenvalue problem, whereas in the following subsections we will illustrate the method of solution 
on two examples. 



A. General considerations 



Plugging Eq. (33) into Eq. (5), we arrive at an ordinary differential equation for (j){p): 

C(j) + E(f)^Q. (35) 

As G{p,t) is an entire function on the complex p-plane [21], (pip) must be analytic in all singular points of differential 
operator C. If the order of this operator is K, this demand yields K "self-generated" boundary conditions for (j){p). 
In view of the equality G{p = l,t) = 1, operator jC vanishes at p = 1, which yields a universal boundary condition: 
(f)(1) = 0. The rest of the K — 1 self-generated boundary conditions are problem-specific, see examples below. 

What is the general structure of differential operator £? For populations that experience extinction, £0 cannot 
include a term proportional to (p, as such a term would correspond to influx of particles into the system, — )• ^, 
and would prevent extinction. In general, C includes first-order derivative terms (corresponding to branching and 
decay processes) and higher-order derivative terms. For extinction scenario A one has '0(0) > 0, see Eq. (32). Let 6o 
denote the rate of decay A — >■ 0, and bm, m = 2, 3, . . . , M, denote the rates of branching reactions A — > mA. One has 
ip^O) = b2 + 2^3 + . . . + (M — 1)&M — bo>0- Rescaling time by V'(O), we see that the (rescaled) coefficient of the term 
n-' (for j = 1,2,. . .) in Eq. (32) must scale as N^~^ to ensure that ni = 0{N). As a result, the (rescaled) coefficient 
of the jth-order derivative term in £ scales as N^~^ , and £ can be written as 

For reaction rates that arc polynomial in n, the functions fj{p) are polynomial in p. Notably, all functions fj{p) vanish 
at p = 1. How does the solution of Eq. (35) look like at A'' 1? As £■ turns out to be exponentially small in A^, the 
simplest approximation for Eq. (35) would be to discard all terms except fi{p)d(p/dp, arriving at a constant solution 
(f>{p) = 1 (according to our choice of normalization). Indeed, as ni = 0{N) ^ 1, the probability to observe n ni 
particles in the metastable state is exponentially small. These probabilities are proportional to low-order derivatives 
oi (p at p — 0, see Eqs. (3) and (33), so (p{p) must indeed be almost constant there. This solution, however, does not 
obey the zero boundary condition at p = 1. The true solution, therefore, must rapidly fall to in a close vicinity of 
p = l, see Fig. 4. The point p = 1 is a singular point of Eq. (35). Actually, when approaching p = 1 from the left, the 
almost constant solution breaks down even earlier: in the vicinity of another point p = Pf < 1 where fi{p) vanishes, 
see the next paragraph. In the vicinity oi p = pf, the first-order derivative term seizes to be dominant, and all terms 
in Eq. (35), including E(p, are comparable. Although (p{p) deviates from a constant value in the vicinity of p = Pf, 
one can still treat this deviation perturbativcly: (p{p) ~ 1 + 5(p{p), where (50 <C 1. When p becomes distinctly larger 
than Pf, (p{p) already varies strongly. Here, the Ecp-teim [which comes from the time derivative of G{p,t)] can again 
be neglected, and so the (nontrivial) solution which is sought in this region is quasi-stationary. The quasi-stationary 
solution can be found in the WKB approximation, as the typical length scale 1/A'', over which (p{p) varies, is much 
smaller here than 1 — pf (a more accurate criterion will appear later). 
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FIG. 4: Shown is a sketch of the eigenfunction (f>{p) of the lowest excited state at A'^ ^ 1 for a typical problem of population 
extinction. is almost constant on the region p < 1 except close to p = 1, where it rapidly goes to zero. 



Why does the root p/ of function exist? After some algebra, function fi{p) can be written as 

M 

^Lb^-P), (37) 

m=0 

where bm = &m/V'(0). The polynomial equation fi{p) = has appeared in the context of n-space description of 
stochastic population extinction [18]. It has been shown in Ref. [18] that this equation has exactly two real roots: 
p ^ 1 and p — pf, where in general < < 1 . 

Now we can summarize the general scheme of solution of the eigenvalue problem for the lowest excited state. One 
has to consider three separate regions: (i) the region to the left, and sufficiently far, from the point p = pf, where one 
can put = 1 up to exponentially small corrections, (ii) in the boundary-layer region \p — Pf\ ^Pf where (f>{p) is 
still very close to 1 and can be sought perturbatively, and (iii) in the quasi-stationary region pf < p < 1, where (f>{p) 
varies strongly, and the WKB approximation can be used. The solutions in the neighboring regions can be matched 
in their joint regions of validity. This procedure holds, at iV 3> 1, for a broad class of systems exhibiting extinction. 

There is a convenient shortcut to this general procedure when the highest-order reaction in the problem is two-body. 
Here the quasi-stationary equation [Eq. (35) with the E(f> term neglected] is always solvable exactly. There is no need 
to apply WKB approximation in such cases, and it suffices to only consider two, rather than three, regions, see the 
next subsection. Finally, regardless of the order of £, it is simpler to deal with u{p) = 4>'{p) rather than with (l){p) 
itself, as this enables one to reduce the order of the ordinary differential equation by one everywhere. 



B. Branching-annihilation-decay 

The first example deals with a population of "particles" which undergoes three stochastic reactions: branching 

A A- 2 A, decay A A and annihilation 2 A A 0. As the state n = is absorbing, the population ultimately 
goes extinct. This example was solved by Kessler and Shnerb [15] via "real-space" WKB approximation, where the 
calculations are done in the space of population size. Here we solve it in the momentum space. Because of the 
presence of the linear decay reaction yl — > 0, this example exhibits a generic transcritical bifurcation as a function 
of the control parameter Rq introduced below, and generalizes simple single-parameter models [16, 19] considered 
earlier. The deterministic rate equation reads 

n = (A — fi)n — af? . (38) 

For A > /i Eq. (38) has, in addition to the trivial fixed point n — 0, also a positive fixed point n\ — (\ — \i)la. When 
starting from any n{t = 0) > 0, the population size flows to the attracting fixed point n = rii, with characteristic 
relaxation time t^. = {\ — /^)~^5 and stays there forever. Rescaling time \t — >■ and introducing rescaled parameters, 
TV = A/cr and — \/ ^jl, the attracting fixed point becomes rii = N{1 — R^^). We demand that ^ 1, and Rq > 1 
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and not too close to 1 (the exact criterion will appear later). When Rq exceeds 1 the deterministic system undergoes 
a transcritical bifurcation. 

To account for intrinsic noise we consider the master equation 

J^Pnit) = ^ [(" + 2)(n + l)Pn+2{t) - n{n - 1)P„(0] + (n - l)P„-i(i) - nP„(i) + + l)fn+i - nP„] , (39) 

where time is rescaled, Xt — )• t. The evolution equation for the probability generating function G{p, t) is 



dt ~ 2N 



l_\ dG 
Ro) dp 



(40) 



At t ^ t.f = (1 — 1/R(f)~^ the metastable probability distribution function, peaked at n ~ ni, sets in, and Eq. (33) 
holds. To determine the QSD and MTE we turn to the Sturm-Liouville problem for the lowest excited eigenmode 
(^(p) and eigenvalue E 



2N 



(i-pV + (p-i) p- 



1 

i?o 



i)' + E(j) = 0. 



(41) 



Here, the self-generated boundary conditions for (t){p) are: (?!)(1) = and 2(1 + i?Q ^)0'(— 1) + E(p{—1) = 0. Because of 
the expected exponential smallncss of E, the latter condition can be safely approximated by 1) ~ 0. 

We now apply the procedure of solution presented in the previous subsection on Eq. (41). Using u{p) = (f)'{p), the 
exact solution of the quasi-stationary equation [Eq. (41) without the E(l) term], 



1 

2N 



(i-pV + (p-i)(p 



1 



u = 0, 



can be written as 



Here 



i(p) = (7e-^^(f). 



S(p) = 2 



1-P+ 1+ 



Rq 



In 



l+p 



(42) 



(43) 



(44) 



To determine the arbitrary constant C we need a boundary condition for u{p) at p = 1. It follows from Eq. (33) that, 



(45) 



On the other hand, by virtue of Eq. (2), the left hand side of Eq. (45) is equal to n{t) which behaves as n\ ex.p{—Et), 
see e.g. Ref. [16]. As a result, u{l) ~ —ni and, by using Eq. (43), we obtain C = -N{1 - i?^ i). Therefore, 



.(p) = -n(i- i-) e-^«(^) 



(46) 



with S{p) from Eq. (44). This yields the solution we looked for: = u{s)ds, which satisfies the boundary condition 
(^(1) = 0. One can check now that neglecting the Ecf) term in Eq. (41) demands pRo — 1 N~^/^. 

Although there is no need in the WKB approximation in this case of a two-body reaction, it is still instructive to 
re-derive Eq. (46) by using the WKB approximation for 0(p). To this end we consider the quasi-stationary version of 
Eq. (41), 



(l-p2)</." + (p-l)(p 



2Ar 



Ro 



0, 



(47) 



and make a WKB ansatz (f>{p) = a{p) cxp[— AfS'(p)]. In the leading order in A'' ^ 1 we obtain a stationary Hamilton- 
Jacobi equation H\p, —S'{p)] = with effective Hamiltonian [25] 



'-rTo 



(1+P)(? 



g(p-i)- 



(48) 
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FIG. 5: Branching-annihilation-decay. Shown are zero-energy lines of Hamiltonian (48) on the (p, q) phase plane. The thick 



solid line corresponds to the instanton q - 
to So from Eq. (55). 



-S'{p) (44). Here qi = m/N — 1 — Rg , and the area of the shaded region is equal 



Here, as in Sec. Ill, q{p) = —S'{p) is the reaction coordinate conjugate to the momentum p. There are two trivial 
zero energy orbits of this Hamiltonian: the deterministic orbit p — I and the "extinction orbit" q = 0. The action 
along the extinction orbit is zero: S{p) = 0, so the corresponding WKB mode can be called "slow". There is also a 
nontrivial zero-energy orbit q{p) = 2{p — i?,^^)/(l -\- p\ It includes a heteroclinic orbit exiting, at t = — oo, the fixed 
point (p = 1, g = gi = ni/N) and entering, at t = oo, the fixed point {p = Rq^ ,q = 0) of the phase plane (p, g), see 
Fig. 5. This orbit is the "extinction instanton" [22, 25]. It describes the most probable path of the system from the 
long-lived metastable state to extinction. Integrating along this orbit and choosing S{p = 1) = 0, we recover Eq. (44). 
This solution can be called the "fast" WKB mode. 

In the subleading order of the WKB approximation one obtains a{p) = (l — i?o ^) (1 +p)/ [2 {p — Rq^)\ for the 
fast, and a{p) = const for the slow WKB modes. The general WKB solution is a superposition of the two modes, 



1 



[l~R^'){l+p) ^ 
2(p-i?„-i) 



-NS{p) 



(49) 



with S{p) from Eq. (44). Here we have already imposed the boundary condition = and normalization condition 
0(0) ~ 1. The p-derivative of (j){p) from Eq. (49) yields, in the leading order, Eq. (46). As it is clear from Eq. (49), 
the WKB solution breaks down in a vicinity of the point p = -RcT^, where the slow and fast WKB modes become 
strongly coupled. Here the quasi-stationarity does not hold. 

We now proceed, therefore, to the non-quasi-stationary region —1 < p ^ Pf (a more restrictive condition will appear 
a posteriori) . It is easier to deal with it in terms of u(p), rather than (/'(p). Here we can treat the E(j) term in Eq. (41) 
perturbatively: 0(p) = 1 -f 5(j){p), where (50 <C 1 [16, 19]. As a result, Eq. (41) becomes an inhomogeneous first-order 
equation for u{p) = 5(j)'{p): 



{p-l)\p-^Ju = -E, 



(50) 



which can be solved by variation of parameter. For two-body reactions the corresponding homogeneous equation, 
which coincides with the quasi-stationary equation (42), is exactly solvable. As a result, one can solve Eq. (50) in the 
entire non-quasi-stationary region which includes both p < pf and \p — pf \ Pf- The solution is 



u{p) 



./-I 1 - 



ds . 



(51) 



where the arbitrary constant is chosen so as to obey the boundary condition m(— 1) — 0. Note, that the integrand in 
Eq. (51) is regular at s = —1, so the perturbative solution is well-behaved. Solution (51) remains valid as long as (f) 
is close to 1. As one can check, this holds for 1 — p ^ N~^^^, cf. Refs. [16, 19]. The perturbative solution (51) can 
be matched with the quasi-stationary solution (46), e.g. at iV"^/^ <C pRo — 1^1 [26]. 
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Solution (51) simplifies in the "left region" p < pf, not too close to pf. By Taylor-expanding the integrand in 
Eq. (51) (which is a monotone increasing ftinction of p for p < Pf) in the vicinity oi s = p, we obtain 



E 



(p_l)(p_i?-i) 



(52) 



This result (which holds in the region 1 — pRo ^ 7V~^/^) has a simple meaning: here the first-derivative term in 
Eq. (50) is negligible. To neglect this term in Eq. (50) [or the term proportional to (^"{p) in Eq. (41)] is the same as to 
disregard the two-body reaction 2 A — > compared with the one-body reactions of branching and decay. This is indeed 
a legitimate approximation at small n [15, 18]. Note that, not too close to p = pf, u{py'^^* is exponentially small in 
A'', so that ~ 1 up to an exponentially small correction. Putting = 1 in the left region, however, would be too a 
crude approximation, as it would only give a trivial left tail of the QSD: tti = tt2 = ■ ■ . = [27]. Correspondingly, the 
solution in the left region cannot be obtained from the WKB approximation. 

We can now find the eigenvalue E by matching the quasi-stationary solution (46) and the perturbative non-quasi- 
stationary solution (51) in their joint validity region N~^/'^ <^ pRo — 1^1 [26]. For pRo — 1 ^ A^~^/^, the integral 



in Eq. (51) can be evaluated by the saddle point approximation. The saddle point is at p : 
is 



Pf 



u{p) 



2EV^RI^'^ 



1 + 



VW+i{Ro - 1) 

Matching this result with the quasi-stationary solution (46), we find 



Rn , and the result 



(53) 



E: 



N{Ro + l){Ro-l)\ _^s. 



An 



r: 



5/2 



(54) 



where 



So = 2 



l-ln2 



ln2 



Rq 



1 

Ro 



In 1 + 



Rf] 



(55) 



The MTE, in physical units, is Tex = with E from Eq. (54), in agreement with Ref. [15]. As i?o — >■ oo the 

decay reaction becomes irrelevant, and one recovers the result for the branching-annihilation model [15, 16, 28]. 
When Ro — 1 <C 1, the system is close to the transcritical bifurcation of the deterministic rate equation. Here the 
Fokker-Planck approximation to the master equation is applicable [13, 15]. The corresponding asymptote of Eq. (54), 



E = 



(i?o-l)'e-^(^''-i)' 



is valid when Rq — 1 ^ N~^/'^, so that E is still exponentially small in A''. 

Having found E, we have a complete solution for u{p), given by Eqs. (46) and (51). Now one can find the QSD by 
using Eq. (3) for n ~ 0(1) and Eq. (4) for n 1. The results coincide with those obtained by Kcsslcr and Shnerb [15] 
by the "real-space" WKB approximation, so we will not present them here. The large-n tail of the QSD decays faster 
than exponentially, thus justifying our a priori assumption that (j){p) is an entire function in the complex p-plane. 
Shown in Fig. 6 is a comparison between the analytical and numerical solutions for dpG c^i —u{p)e~^* at a time 

< < < 1/-E, when dpG ~ -u{p). 



C. Branching and triple annihilation 



Here we again consider a metastable population on the way to extinction, but now a three-body reaction is present. 
Our model system includes two reactions: the branching A^2A and the triple annihilation 3^ % 0. The deterministic 
rate equation, 

^ = An-^#, (56) 

has two relevant fixed points: the repelling point n = and the attracting point ni = (2A//i)^/^ = iV ^ 1. According 
to Eq. (56), the system size approaches n = ni after the relaxation time tr = A^^, and stays there forever. Contrary 



14 




FIG. 6: Branching-annihilation-decay. Shown is the p-derivative dpG of the probability generating function G at tr ^ f <C 1/E 
for A'^ = 10'^ and Ro = 1.5. The solid line is the absolute value of the perturbative solution (51); the dashed line is the 
absolute value of the quasi-stationary solution (46). In the joint region of their validity the two lines are indistinguishable. The 
crosses indicate the values obtained by a numerical solution of Eq. (40) for no = 20 particles at f = and boundary conditions 
G(l, t) = 1 and dpG{-l, t) = Q. 



to this prediction, fluctuations drive the population to extinction. Upon rescaling time t — > At, the master equation 
reads 

= (n - l)P„_i - nP„ + ^ [(n + 3)(n + 2){n + l)P„+3 - n{n - l){n - 2)P„] , 

(57) 

whereas the evolution equation for Gij), t) is 

At t t,., Eq. (33) holds, and the ordinary differential equation for the lowest excited eigenfunction <j){p) is 

^ (1 - p^)(^"' + Pip - + ^ . (59) 



3A^2 



This equation has three singular points in the complex p-plane. These are the roots of 1 — p'^: one real, pi = 1, and 
two complex, p2 = e^'^*/'^ and = e^^'/'^. Since (/)(p) must be analytical in all these points, (/)(p) must satisfy three 
conditions: 

p,{p,-l)q^'{p,) + Ecb{pi)^0 1,2,3. (60) 

Here the p-derivative is in the complex plane. For i = 1 Eq. (60) yields <f>{p = 1) = 0. As turns out to be 
exponentially small in TV ^ 1, we can neglect small terms proportional to E in the conditions for i — 2 and 3 and 
obtain 0'(p = e^'''/^) ~ and (j)'{p = e'''^*/^) ~ 0. 

In the quasi-stationary region (the exact location of which will be determined later) Eq. (59) becomes 

^ (l-p3)0"'-Fp(p-l)0' = O. (61) 



3iV' 



This equation is of second order for u(p) = (j)' {p), but it is not exactly solvable in terms of known special functions, 
and this is a typical situation for three-body, four-body, . . ., reactions. The presence of the large parameter N ^ 1 
justifies the WKB ansatz 0(p) = a(p)e~^'^'-^-'. It yields, in the leading order of iV ^ 1, a stationary Hamilton- Jacobi 
equation H[p, — S"(p)] = with Hamiltonian 

H{p,,)=\p- ^'^'y'^'' ],{p-l). (62) 
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FIG. 7: Branching and triple annihilation. Shown are zero-energy lines of Hamiltonian (62) on the {p, q) phase plane. The 
thick solid line is the instanton q — ~S'{p) = i'ip), given by Eq. (63). Here qi — 1, and the area of the shaded region is equal 
to So from Eq. (75). The dashed line denotes a non- physical orbit. 



Here again, in addition to the trivial zero-energy lines q = and p = 1, one obtains an instanton orbit 

q^i,{p)^(-J^Y' (63) 



which connects the fixed points (1, qi = ni/N = 1) and (0, 0) in the {p, q) plane, see Fig. 7. The instanton corresponds 
to the fast-mode WKB solution, whereas the orbit q — corresponds to the slow-mode WKB solution, similarly to 
the previous example. 

Again, it is simpler to do the actual calculations for u{p) — (t)'{p), rather than for 4){p). Using the WKB ansatz 
u{p) = b{p)e~^^^^^ in the quasi-stationary equation 

^ {l+p + p^y -pu^O, (64) 



37V' 
we obtain 



jl+p + p^) 
3A^2 



[N^iS' fb - 2NS'b' - NS"b] -pb = 0, (65) 



where we have neglected the sub-subleading term proportional to b" /N'^. In the leading order we obtain S'{p) — —ip{p) 
[the solution with S'{p) = ip{p) is non-physical and must be discarded]. The arbitrary constant can be fixed by putting 
5(1) = 0, and we obtain 

S{p) ^ - J\{x)dx, (66) 

with ■ijj{x) from Eq. (63). This result can be expressed via elliptic integrals, but we will not need these cumbersome 
formulas. 

In the subleading order Eqs. (65) and (66) yield a first-order ordinary differential equation for b{p) whose general 
solution is 

C(l+p + p^)V^ 

^(P^ = ^173 ■ (67) 

We demand u(l) ~ — rti = —N [see Eq. (45)] and obtain the quasi-stationary WKB solution for u{p): 

u-^-^ip) = -^^^^^^0^e~-s^^^ (68) 
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with S{p) from Eq. (66). Now one can check that asymptote (68) is vahd when p ':$> N~'^/'^] otherwise it is not 
justified to neglect the term h" [p) in Eq. (65). Again, the quasi-stationarity and the WKB approximation break down 
in a vicinity of the point where the fast and slow WKB modes are strongly coupled. In this example this point is at 
p = 0, whereas in the previous example it was at p = p/ ^0. That the WKB breaks down here at p = is a special, 
non-generic situation resulting from the absence of the linear decay process A — >■ from the set of reactions A — > 2A 
and 3^ — > we are dealing with. 

To remedy the divergence of the WKB solution aip = 0, one need to account for a deviation from quasi-stationarity. 
The corresponding non-quasi-stationary solution of Eq. (59) is perturbative in E, as in the previous example, so the 
equation we need to solve is 

^ {l-p^)u" +p{p-l)u = -E. (69) 



The corresponding homogeneous equation, Eq. (61), is not solvable in known special functions. Therefore, we will 
solve Eq. (69) approximately in two separate regions and match the solutions in their joint region of validity. 

The first region, which we call "left", is p < (and not too close to zero, see below). Here we can neglect the 
u"-term in Eq. (69) and obtain 

"'■"<^' = Kl^- <™> 

This asymptote, valid when ~p ^ 7V~^/^, corresponds to neglecting the high-order reaction iA — >■ at small 
population sizes. As in the previous example, u''^'^*(p) is exponentially small. By choosing an exponentially small 
solution for u{p) in the left region, we effectively discarded two other linearly independent solutions of Eq. (61) which 
are singular at p = e^'^*/'^ and p = e^'^^l'^ . As pf = here, one can actually put w = in the left region and still 
accurately determine the QSD [27]. 

The second region is the boundary layer |p| <C 1, where Eq. (69) becomes 

^ u"-pu = -E, (71) 



3Af2 



The general solution of this equation is 



u'\p) 



rP rP 
Ci+a^nE / Bi{as)ds Ai{ap) + C2 — oP"kE i Ai{as)d. 
Jo \ V Jo 



Bi{ap) , (72) 



where Aiijj) and Bi{y) arc the Airy functions of the first and second kind, respectively [23], and a = {iN'^Y/^ . 

Now we can find the unknown constants ci and C2 (assuming for a moment that E is known) by matching the 
asymptotes (70) and (72) in their common region N~'^/^ <i; — p <C 1. As w'^-''*(p) is exponentially small at AT^Ip]^ ^ 1, 
the boundary layer solution u''\p) from Eq. (72) must also be exponentially small there. Evaluating the integrals in 
Eq. (72) at p = — oo and using the identities Bi{s)ds = and Ai{s)ds = 2/3, we arrive at 

2ttEN^/^ 

c,c,0, C2C ^j^. (73) 

Now we can find the extinction rate E by matching the asymptotes of u^^^{p) and ul'^p) in their common region 
jY-2/3 ^p^i_ The p < 1 asymptote of the WKB solution (68) is 

" - (3p)l/4^ ^ ' ^^^> 

where 

3:. 



So = J^ ^TTTT^j ^- = 0-836367..., (75) 
is the shaded area in Fig. 7. Let us obtain the p » A^-^/s asymptote of u'''(p) (72). First, for z » 1 [23] 
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FIG. 8: Branching and triple annihilation. Shown is a comparison between the extinction rate (78) (solid line) and the extinction 
rate — [ln(l — Po^"™(t))]/f (crosses) found from a numerical solution of the master equation (57) at different N. The inset shows 
the ratio of the two rates. 



Now we need to evaluate the integrals in Eq. (72). As we are interested in the region of N^p'^ ^ 1, the integral of 
Ai{as) can be evaluated by putting p = oo and using the saddle point approximation, arriving at 



Ai 



2x1/3, 



ds 



34/3^2/3 



The main contribution to the integral of Bi{as) at N'^p'^ ^ 1 comes from a vicinity oi s = p, where Bi{as) is 
exponentially large, see Eq. (76). Expanding the exponent in a Taylor series around s ~ p, we obtain in the leading 
order 



p 



Bi 



(37V2)l/3^ 



ds 



,(2/V3)JVp3/=* 



37/127^1/2 jV7/6p3/4 ■ 

Now one can see from Eq. (73) that the main contribution to u^\p) (72) comes from the Bi{ap) term, and we obtain 
Matching Eqs. (74) and (77), we obtain 



i? = y^e-^^«. (78) 

The MTE in physical units is given by t^x = {XE)~^, which is exponentially large in N, as expected. A comparison 
between the analytical result for the extinction rate (78) and a numerical result, obtained by solving (a truncated 
version of) master equation (57), is shown in Fig. 8. For ^ 1 the agreement is excellent. 
Now let us calculate the QSD. Combining Eq. (4) with Eqs. (33) and (34), we obtain 

n,,>^ = --^<f^dp. (79) 

For n 1 we can use the WKB asymptote (68): 



27rm J p"- 2TTni J (3]3)i/4 



with given by Eq. (63). As A^ 3> 1 and n ^ 1, this integral can be evaluated via saddle point approximation 
[24]. Let us denote f{p) — A^ ^{x)dx — n\np. The saddle point equation f'{p*) — reduces to a cubic equation 
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which has one and only one real root = (n/N). As f"{p*) > 0, we must choose a contour in the complex p-plane 
which goes through this root perpendicularly to the real axis. The Gaussian integration yields 



\ +p,+ piy/"^ exp [N J[' ^{x)dx] 
n^2^/"(p,)(3p,)i/4 P2 



(82) 



we omit a cumbersome expression for f"(p)- Note, that for n ^ 1, the saddle point is always obtained in the 
region where u^^^{p) is valid, see below. Let us calculate the 1 <^ n N and n ^ N asymptotes of Eq. (82) 
with an exponential accuracy, ln7r„ ~ f{p*)- For n <^ N the saddle point, given by Eq. (81), is obtained at 
= [n/ (VSA^)]^/'^ <C 1. Here it suffices, in the leading order in n/N, to put = in the upper bound of the integral 
in Eq. (66). Then the integral yields 5*0 from Eq. (75). For n':$> N we obtain = [n/{V3N)]'^ ^ 1. Here a dominant 
contribution to the integral in Eq. (66) comes from the region of ^ 1 which enables one to simplify the integrand. 
The resulting asymptotes are 



ln7r„ ~ N 
IniTn ^ N 



2n N n\ 
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n 
N 



n 
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0{l) 



, n:^N. 



(83) 



Notice that each of these tails of the QSD are non-Gaussian. The n ^ N tail decays faster than exponentially, thus 
justifying a posteriori our assumption that (/)(p) is an entire function in the complex p-plane. 

At \n — N\ -C iV the saddle point, given by Eq. (81), is obtained at p* = 1 + {n — N)/N, and we arrive at a Gaussian 
asymptote 



1 



V2^ 



-{n-NY /{2N). 



(84) 



the preexponent is fixed by normalization. Equation (84) holds for |n — iV| ^ N^^^; this condition is tighter than 
\n — N\ <^ N. Note, that the Gaussian asymptote of the QSD can also be found by directly calculating the mean and 
variance of the distribution. These (and other higher cumulants of the distribution) can be found by using derivatives 

of G{p,t) with respect to p at p = 1, see Eq. (2). Indeed, from Eqs. (2) and (6), the mean of the QSD (at times 



< i < Tecc) is given by n = dpG\p=i ~ -u(p = 1) 
In its turn the variance in the leading order is 



N, where here we have used 



,WKB 



(p) given by Eq. (68). 



'Pnit) - 5]nP„(t) = [dppG + dpG-{dpGf] ~ -u'{l) - u{l) - [u{l)f ~ AT, 



n=0 



recovering the Gaussian asymptote (84). 

At n = 0{1) the QSD can be evaluated directly from 



TTn = — 



1 rf"-iu(p) 



n! dp 



,71—1 



n> 1. 



p=0 



Here one should use the boundary-layer solution around p = 0, given by Eqs. (72) and (73). This yields 



r(i/3)£Ar2/3 

"3V3 



772 = 



3i/6r(i/3) 



and TTs 



(85) 



To calculate other n = 0(1) terms, one can use a recursion relation obtainable from the master equation (57) with 
P„ = 0. Indeed, at n <C one can neglect the terms n(n — l)(n — 2)P„ and (n — l)P„_i compared with the terms 
(n + 3)(n + 2)(n + l)P„+3 and nP„, respectively, and arrive at the following relation: 

= 7 — -^;t7 — -^TTT — rrrrT^n ■ (86) 



(n + 3)(n + 2)(n+l) 



Note, that the small-n (86) and the WKB (82) segments of the QSD have a joint region of validity at 1 <C n ^ A^. 

A comparison between WKB result (82) and a numerical solution of (a truncated version of) master equation (57) 
is shown in Fig. 9. The inset compares the n <^ N analytical asymptote [see Eqs. (85) and (86)] with numerical 
results. Excellent agreement is observed in both cases. It can be also seen that the Gaussian approximation (84) 
strongly overestimates the QSD in the low-n region, and underestimates it in the high-n region. 
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FIG. 9: Branching and triple annihilation. Shown is the natural logarithm of the QSD versus n for A'^ = 20. The dashed line 
is WKB solution (82), the dash-dotted line is the Gaussian approximation (84), and the solid line is the numerical solution of 
the (truncated) master equation (57). Inset: the n <^ N asymptote of the QSD obtained analytically [Eqs. (85) and (86)] (x's) 
and numerically (fat dots). 

V. DISCUSSION 

The p-space representation renders a unique perspective to theory of large fluctuations in populations undergoing 
Markovian stochastic gain-loss processes. The stationary distribution of the population size is encoded in the ground- 
state eigenfunction of a Sturm-Liouville (spectral) problem for the probability generating function. In the case 
of a long-lived metastable population on the way to extinction, the MTE and the quasi-stationary distribution of 
population size are encoded in the eigenfunction of the lowest excited state. The uniqueness of solution in these 
problems is guaranteed by the condition that the probability generating function is an entire function on the whole 
complex p-plane except at infinity. As this work has demonstrated (see also Refs. [16, 19, 20, 22, 25, 29]), the p-space 
representation in conjunction with the WKB approximation and other perturbation tools employing a large parameter 

^ 1 (the mean population size in the stationary or metastable state) yields accurate results for extreme statistics 
in a broad class of problems of stochastic population dynamics. Such an accuracy is usually impossible to attain via 
the van Kampen system size expansion which approximates the exact master equation by a Fokker-Planck equation. 

How does the p-space approach compare with the "real" space WKB method of Refs. [5, 6, 15, 17, 18] when the 
stationary or metastable population size is large, iV ^ 1? One advantage of the p-space representation is that, for 
two-body reactions, there is no need in the WKB approximation, as the quasi-stationary equation in this case is 
always solvable exactly. Another advantage appears when the WKB solution for G{p, t) is valid for every p ^ 0, as 
occurs in the molecular hydrogen production problem. Sec. III. In such cases one directly finds the entire probability 
distribution function, including the region of small n = 0(1). In the real-space approach a separate (non-WKB) 
treatment of the n = C(l) region, and a matching with the WKB-solution valid at n ^ 1 would be needed [18]. 

Still, from our experience, every problem which includes the large parameter A^ 3> 1, and can be solved in the p- 
space, can be also solved in the "real" space. Furthermore, for populations exhibiting escape to infinity [5], escape to 
another metastable state [17], or Scenario B of extinction [18], the p-space representation meets significant difficulties. 
One difficulty is that one should account for a constant-current WKB solution in these cases [5, 17, 18]. The constant- 
current solution comes from the deterministic line p = 1 of the phase plane of the underlying classical Hamiltonian. In 
the ^-representation this line is vertical, as in Figs. 1, 5 and 7, and so the constant-current solution cannot be easily 
accounted for. In addition, it is unclear how to deal with the region of non- uniqueness oi q = q{p) which is inherent, 
in the p-representation, in these cases. There are two WKB solutions in this region, one of them exponentially small 
compared with the other. The real-space approach avoids these difficulties, and the solution in these cases can be 
worked out in a straightforward manner [5, 17, 18]. 

An important advantage of the p-space representation stems from the fact that the evolution equation for G(p, t) 
is exactly equivalent to the original master equation. Therefore, the p-space approach is especially valuable for exact 
analysis, as illustrated by the example of molecular hydrogen production, see Ref. [8] and Section III. 

Finally, generalization of the p-representation to interacting multi-species populations is quite straightforward, see 
Ref. [30] . The resulting multi-dimensional evolution equation for the probability generating function can be analyzed in 
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WKB approximation. As of present only the leading-order WKB-approximation for population extinction is available, 
and this is regardless of whether one uses the p- or n-space approach. In the leading WKB order the problem again 
reduces to finding a nontrivial zero-energy trajectory of the corresponding classical Hamiltonian, and the action along 
this special trajectory. This problem can be solved numerically. If additional small parameters are present, the 
problem may become solvable analytically, again in both p- and n-spaces [30-32] . 
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